/* 3D Navier-Stokes Solver <-> (Level Set <-> Fast March) -> Marching Cubes */
#include <iostream>
#include "typen.h"
#include "matrix.h"
#include "navier.h"
#include "convective.h"
#include <math.h>
#include <errno.h>
#include <stdlib.h>
#include <cstring>
#include "fastMarch.h"

using namespace std;
int  neighbors_table[64][7] ={{0,-1,-1,-1,-1,-1,-1},{1,0,-1,-1,-1,-1,-1},{1,1,-1,-1,-1,-1,-1},{2,1,0,-1,-1,-1,-1},
                               {1,2,-1,-1,-1,-1,-1},{2,2,0,-1,-1,-1,-1,},{2,2,1,-1,-1,-1,-1},{3,2,1,0,-1,-1,-1},
                               {1,3,-1,-1,-1,-1,-1},{2,3,0,-1,-1,-1,-1},{2,3,1,-1,-1,-1,-1},{3,3,1,0,-1,-1,-1},
                               {2,3,2,-1,-1,-1,-1},{3,3,2,0,-1,-1,-1},{3,3,2,1,-1,-1,-1},{4,3,2,1,0,-1,-1},
			       {1,4,-1,-1,-1,-1,-1},{2,4,0,-1,-1,-1,-1},{2,4,1,-1,-1,-1,-1},{3,4,1,0,-1,-1,-1},
                               {2,4,2,-1,-1,-1,-1},{3,4,2,0,-1,-1,-1},{3,4,2,1,-1,-1,-1},{4,4,2,1,0,-1,-1},
                               {2,4,3,-1,-1,-1,-1},{3,4,3,0,-1,-1,-1},{3,4,3,1,-1,-1,-1},{4,4,3,1,0,-1,-1},
                               {3,4,3,2,-1,-1,-1},{4,4,3,2,0,-1,-1},{4,4,3,2,1,-1,-1},{5,4,3,2,1,0,-1},
			       {1,5,-1,-1,-1,-1,-1},{2,5,0,-1,-1,-1,-1},{2,5,1,-1,-1,-1,-1},{3,5,1,0,-1,-1,-1},
                               {2,5,2,-1,-1,-1,-1},{3,5,2,0,-1,-1,-1},{3,5,2,1,-1,-1,-1},{4,5,2,1,0,-1,-1},
                               {2,5,3,-1,-1,-1,-1},{3,5,3,0,-1,-1,-1},{3,5,3,1,-1,-1,-1},{4,5,3,1,0,-1,-1},
                               {3,5,3,2,-1,-1,-1},{4,5,3,2,0,-1,-1},{4,5,3,2,1,-1,-1},{5,5,3,2,1,0,-1},
			       {2,5,4,-1,-1,-1,-1},{3,5,4,0,-1,-1,-1},{3,5,4,1,-1,-1,-1},{4,5,4,1,0,-1,-1},
                               {3,5,4,2,-1,-1,-1},{4,5,4,2,0,-1,-1},{4,5,4,2,1,-1,-1},{5,5,4,2,1,0,-1},
                               {3,5,4,3,-1,-1,-1},{4,5,4,3,0,-1,-1},{4,5,4,3,1,-1,-1},{5,5,4,3,1,0,-1},
                               {4,5,4,3,2,-1,-1},{5,5,4,3,2,0,-1},{5,5,4,3,2,1,-1},{6,5,4,3,2,1,0}};

bool update_pressure = false;
void Navier::Navier_Solver_Initialization() 
//void Navier::DefinirCenaInicial(Matrix <flagtype> flag) 
{ 
   //cout<<"\n Entrei na Definir Cena Inicial \;
        cout<<":Navier_Solver_Initialization \n ";
     
       num_iter_Navier=0;
       dt = 0;
       deltaUVW = size_vectors;
       for(int it=0;it<=5;it++)
	  container_type[it] = it + 1;
       P = new double[size_vectors];
       fill(P,P+size_vectors, 0); 
         
       RHS= new double[size_vectors];
       fill(RHS,RHS+size_vectors, 0);
         
       BiCG_r0 = new double[size_vectors];
       fill(BiCG_r0, BiCG_r0+size_vectors, 0);
       output_area[0] = BiCG_r0;
       
       BiCG_rj = new double[size_vectors];
       fill(BiCG_rj, BiCG_rj+size_vectors, 0);
       output_area[1] = BiCG_rj;
       
       BiCG_pj = new double[size_vectors];
       fill(BiCG_pj, BiCG_pj+size_vectors, 0);
       output_area[2] = BiCG_pj;
       
       BiCG_vj = new double[size_vectors];
       fill(BiCG_vj, BiCG_vj+size_vectors, 0);
       output_area[3] = BiCG_vj;
       
       BiCG_sj = new double[size_vectors];
       fill(BiCG_sj, BiCG_sj+size_vectors, 0);
       output_area[4] = BiCG_sj;
       
       BiCG_tj = new double[size_vectors];
       fill(BiCG_tj, BiCG_tj+size_vectors, 0);
       output_area[5] = BiCG_tj;
       
       MatBuf = new double[size_vectors];
       fill( MatBuf,  MatBuf+size_vectors, 0);
       cout<<"Definir Cena Inicial - All Poisson Matrixes Initialized\n ";
	
      int djx4 = 4 * dj;
      int index = 2*(dk+dj+1);
      for(int k=2;k <= nz+1;k++)
      {
         for(int j=2;j <= ny+1;j++)
	 {
            for(int i=2;i <= nx+1;i++)
	    {
               flag[index] = AIR;
	       index++;
	    }
	    index +=4;
	 }  
	 index += djx4;  
      }       
      int  dj_1 = dj - 1;
      for(int i= 0;i < size_vectors;i +=dj)
      {   flag[i] += 1;flag[i+dj_1] += 2;}
      int dk_dj = dk - dj;
      for(int i=0; i<= nx+3; i++)
      {
          for(int j=i; j < size_vectors; j +=dk)
	  {   flag[j] += 3;flag[j+dk_dj] += 6;}
      }
      int size_vectors_dk = size_vectors - dk;
      for(int i=0; i<= nx+3; i++)
      {
          for(int j=0; j <= ny+3; j ++)
	  {   int k=i + j*dj; flag[k] += 9;flag[k+size_vectors_dk] += 18;}
      }
       cout<<"Definir Cena Inicial - Wall around recipient constructed\n ";
      Init_inflow_and_Label_Slip_Outflow();
      cout<<"Definir Cena Inicial - Initial cell values determined\n ";
      S->DumpInfo(stdout);
      cout<<"Definir Cena Inicial - Parameters shown \n ";
}
void Navier::Init_inflow_and_Label_Slip_Outflow()
{
   std::vector <int> List_Inflow;
   std::vector <int> List_Slip_Outflow;
   int index= dk + dj + 1;
   int djx2 = 2 * dj;
   for(int k=1;k<=nz+2;k++)
   {
      for(int j=1;j<=ny+2;j++)
      {
         for(int i=1;i<=nx+2;i++)
         {
	    if(gridPhi[index] <= 0) //it is inflow
            {
               flag[index] = INFLOW;
	       for(int lt=0; lt<=2;lt++)
	       {  
                   gridUVW[lt][index] = inflow_value[lt]; 
		   gridCellUVW[lt][index] = inflow_value[lt];
		   gridUVW[lt][index+size_vectors] = inflow_value[lt];
	       }
               List_Inflow.push_back(index);
            }
	    else
	    { 
	     	if(ISOUTFLOWIJK(i,j,k)) 
            	{
		   //cout<<"outflow encontrado ="<<index<<"\n ";
		   //int a; cin>>a;
                   flag[index] = OUTFLOW_AIR;
		   List_Slip_Outflow.push_back(index);
            	}
	    	else 
		{
		   if(ISSLIPIJK(i,j,k))
		   //if((i==1) || (j==1) || (k==1) || (i==S->gridp[0]+1) || (k==S->gridp[2]+1)) 
	           {
                      flag[index] = SLIP_AIR;
                      List_Slip_Outflow.push_back(index);
                    }
		 }
              }
	    index ++;
         }
	 index +=2;
      }
      index += djx2; 
   }
   
    int itend = (int)List_Inflow.size();
    for(int itfor = 0; itfor < itend; itfor ++)
    {
    	Initialize_Border_Faces_Lists(List_Inflow[itfor]);
    }
    itend = (int)List_Slip_Outflow.size();
    for(int itfor = 0; itfor < itend; itfor ++)
    {
    	 Label_Slip_Outflow_Cells(List_Slip_Outflow[itfor]);
    }
    List_Border_Face_size_m1 =  List_Border_Face_Air->size()-1;
   for(int jt=0;jt<=2;jt++)
       dt = fmin(dt, h/(fabs(inflow_value[jt])+1E-15));
}

void Navier::Label_Slip_Outflow_Cells(int index)
{ 
      static int d[6] = {-1,1,-dj,dj,-dk,dk};
      static int power[6] = {1,2,4,8,16,32};
      int neighbors_free = 0;
      double power_sum = 0;
      int n1;
     // cout<<"Label_Slip_Outflow_Cells - Starting point\n ";
      for(int n=0; n<6; n++)
      {
         int indexn = index +d[n];
         if(flag[indexn] == AIR)
         {
	    neighbors_free++;
	    power_sum += power[n];
            if(neighbors_free == 1)
            {
                  n1 = n;
            }
               
         }
      }
      if(neighbors_free == 1)
      {
         flag[index] += incrementby4 + n1;
      }
      if(neighbors_free > 1)
      {
         flag[index] += power_sum;
      }
      
     // cout<<"Label_Slip_Outflow_Cells - Ending point\n ";
}

void Navier::Initialize_Border_Faces_Lists(int index)
{ 
            static int jump[6]={1,dj,dk, -1,-dj,-dk};
	    cout<<"Initialize_extended_border_Lists - Starting point"<<index<<"\n ";
            for(int it= 0; it<=2;it++)
	    {
	        int index_neighbor = index+ jump[it];
		if(( flag[index_neighbor]>= AIR) && (flag[index_neighbor] != INFLOW))
		{
		     gridUVW[it][index_neighbor] = inflow_value[it];
		     gridUVW[it][index_neighbor+size_vectors] = inflow_value[it];
		     *Total_inflow += inflow_value[it];
		     List_Border_Face_Air->push_back(index_neighbor);
                     List_Border_Face_Fluid->push_back(index);
		    
		}
		 
	    }
	    for(int it= 3; it<=5;it++)
	    {
	        int index_neighbor = index+ jump[it];
		 
                if(( flag[index_neighbor]>= AIR) && (flag[index_neighbor] != INFLOW))
		{
		    *Total_inflow += inflow_value[it - 3];
		    List_Border_Face_Air->push_back(index_neighbor);
                    List_Border_Face_Fluid->push_back(index); 
		}
		
		
	    }
	    cout<<"Initialize_extended_border_Lists - Ending point\n ";
}
         
            
void Navier::Processing_Border_Air_List()
{
    cout<<"Processing_Border_Air_List - starting point \n";
    static int dlocal[6]={1,-1,dj,-dj,dk,-dk};
    static int face_bit[6] ={1,2,4,8,16,32};
    vector<int>::iterator itair = List_FluidToAir.begin();
    int *antecessor= new int[List_Border_Air->size()];
    fill(antecessor,antecessor+List_Border_Air->size(),0.0);
    int *antecessor_border= new int[List_Border_Air->size()];
    fill(antecessor_border,antecessor_border+List_Border_Air->size(),0.0);
    int *pointer= new int[List_Border_Air->size()];
    int *back_pointer= new int[List_Border_Air->size()];;
    int processed_number = 0;
    int curr_position =0;
    int curr_index;
    int itend = (int)List_Border_Air->size();
    for(int itfor = 0; itfor < itend; itfor ++)
    {
        back_pointer[itfor]=itfor-1;
        pointer[itfor] = itfor +1;
    }
    //cout<<"Processing_Border_Air_List - pointers structure mounted \n";
    pointer[List_Border_Air->size() - 1] = 0;
    back_pointer[0]= List_Border_Air->size()-1;
    
    for(curr_position = 0; curr_position < itend; curr_position++)
    { 
       curr_index = (*List_Border_Air)[curr_position]; 
       if(flag[curr_index]> PROCESSED_BORDER_AIR)
       {
            //cout<<"Processing_Border_Air_List - slip or outflow border cell found \n";
            pointer[back_pointer[curr_position]] = pointer[curr_position];
            back_pointer[pointer[curr_position]] = back_pointer[curr_position];
	    //cout<< "%%%%%%%%%%%%curr_position - antecessor_border[curr_position]%%%%%%%%" << curr_position<<","<< antecessor_border[curr_position] <<"\n";
	    processed_number++;  
	    List_Border_Slip_Outflow_Air.push_back(curr_index);
	    //cout<< "curr_index - flag[index]=" << curr_index<<","<< flag[curr_index] <<"\n"; 
	     
       }
       else
       {
          int sign =-1;
          for(int ifor=0;ifor < 6; ifor ++)
          {   
             int curr_index_temp = curr_index + dlocal[ifor];
             if(flag[curr_index_temp] >= BORDER_AIR)
             {  
	        int jfor = ifor/2; 
                double gridCellUVW_temp = gridCellUVW[jfor][curr_index_temp];
                if((gridCellUVW_temp + gridCellUVW[jfor][curr_index]) * sign > 0)
		{
                      antecessor[curr_position] += face_bit[ifor];
		      if(flag[curr_index_temp] == BORDER_AIR)
		      {
		         antecessor_border[curr_position] += face_bit[ifor];
		      }
	        }  
	     }
	     sign = -sign;
          }
	  //cout<<"Processing_Border_Air_List - a new cell processed \n";
          if(antecessor_border[curr_position] == 0)
          {
             Ordered_New_Fluid_Cell_Initializer(curr_index,neighbors_table[antecessor_border[curr_position]][0]);
             flag[curr_index]= PROCESSED_BORDER_AIR;
	     //cout<<"curr_position ="<<curr_position<<"\n";
	     //cout<<"pointer[curr_position] ="<<pointer[curr_position]<<"\n";
	     //cout<<"back_pointer[curr_position] ="<<back_pointer[curr_position]<<"\n";
             pointer[back_pointer[curr_position]] = pointer[curr_position];
             back_pointer[pointer[curr_position]] = back_pointer[curr_position];
             processed_number++;
	     //cout<< "curr_position - antecessor_border[curr_position]" << curr_position<<","<< antecessor_border[curr_position] <<"\n";
          } 
	  else
	  {
	     //cout<< "*****curr_position - antecessor_border[curr_position]******" << curr_position<<","<< antecessor_border[curr_position] <<"\n";
	  }
      }
    }
      //cout<<"Processing_Border_Air_List - First Turn Completed \n";
      int itendm1 =itend -1;
      curr_position = itendm1;
      //cout<< "processed_number,itendm1 ="<< processed_number<<","<<itendm1 << "\n";
      //exit(0);
      while(processed_number < itendm1)
      {  
      //cout<< "processed_number,itendm1 ="<< processed_number<<","<<itendm1 << "\n";
           int *ant_number_pointer, *ant_pointer;
	   int ant_processed_number = processed_number;
	   do
	   {  
	        curr_position = pointer[curr_position];
		//cout<< "curr_position - antecessor_border[curr_position] \n";
		//cout<< "processed_number,itendm1 ="<< processed_number<<","<<itendm1 << "\n";
		//cout<< curr_position<<","<< antecessor_border[curr_position] << "\n";
	        curr_index = (*List_Border_Air)[curr_position];
		//cout<<"curr_index,flag[curr_index]="<<curr_index<<","<<flag[curr_index]<<"\n";
		//int x,y,z;
		//GIJK(curr_index,x,y,z);
                //cout<<"x,y,z="<<x<<","<<y<<","<<z<<"\n";
		////cout<<"antecessor_border[curr_position]="<<antecessor_border[curr_position]<<"\n";
		ant_number_pointer = &(neighbors_table[antecessor_border[curr_position]][0]);
		ant_pointer = ant_number_pointer + 1; 
		//cout<<"Processing_Border_Air_List - Pointers Border Table set: ant_number_pointer,ant_pointer  \n"; 
		//cout<< ant_number_pointer<<","<< ant_pointer << "\n";  
		//cout<< "*ant_number_pointer ="<< *ant_number_pointer << "\n";
		
	        for(int ifor=1;ifor <= *ant_number_pointer; ifor++)
	        {
	            int neighbor_temp = *ant_pointer;
		    //cout<< "neighbor_temp ="<< neighbor_temp << "\n";
                    int curr_index_temp = curr_index + dlocal[neighbor_temp];
                    if(flag[curr_index_temp] != BORDER_AIR)
                    {   
                         antecessor_border[curr_position] -= face_bit[neighbor_temp];
	            }
	            ant_pointer ++;
		    //cout<< *ant_pointer << "\n";
		    
	        }
		//cout<<"Processing_Border_Air_List - Other Turns new cell processed \n";
                if(antecessor_border[curr_position] == 0)
                { 
                    Ordered_New_Fluid_Cell_Initializer(curr_index,*ant_number_pointer);
                    flag[curr_index]= PROCESSED_BORDER_AIR;
	            //cout<<"curr_position ="<<curr_position<<"\n";
	            //cout<<"pointer[curr_position] ="<<pointer[curr_position]<<"\n";
	            //cout<<"back_pointer[curr_position] ="<<back_pointer[curr_position]<<"\n";
                   pointer[back_pointer[curr_position]] = pointer[curr_position];
                   back_pointer[pointer[curr_position]] = back_pointer[curr_position];
                   processed_number++;
                } 
	    }
	    while(curr_position < pointer[curr_position]);  
	    //cout<<"Processing_Border_Air_List - A turn more processed \n";   
            if(processed_number == ant_processed_number)
	    {
                int ant_number, ant_number_min, curr_position_min;
		
		 curr_position = pointer[curr_position];
		 curr_position_min = curr_position;
		 ant_number_min = neighbors_table[antecessor_border[curr_position]][0];
		
		while(curr_position < pointer[curr_position])
	        {  
	             curr_position = pointer[curr_position];
		     ant_number = neighbors_table[antecessor_border[curr_position]][0];
		     if(ant_number < ant_number_min)
		     {
		        ant_number_min = ant_number;
			curr_position_min = curr_position;
		      }
		}
	       // cout<<"Processing_Border_Air_List - minimum antecessors cell found \n";   
		int curr_index_min = (*List_Border_Air)[curr_position_min];
		//cout<<"curr_index_min,flag[curr_index_min]="<<curr_index_min<<","<<flag[curr_index_min]<<"\n";	
                Ordered_New_Fluid_Cell_Initializer(curr_index_min,ant_number_min);
                flag[curr_index_min]= PROCESSED_BORDER_AIR;
	        //cout<<"curr_position ="<<curr_position<<"\n";
	       // cout<<"pointer[curr_position] ="<<pointer[curr_position]<<"\n";
	        //cout<<"back_pointer[curr_position] ="<<back_pointer[curr_position]<<"\n";
                pointer[back_pointer[curr_position_min]] = pointer[curr_position_min];
                back_pointer[pointer[curr_position_min]] = back_pointer[curr_position_min];
                processed_number++;
            }
	    // cout<<"Processing_Border_Air_List - ending while inside\n";       
      }
      //cout<<"Processing_Border_Air_List - ending while outside \n";
      Ordered_New_Fluid_Cell_Initializer(curr_index,neighbors_table[antecessor_border[curr_position]][0]);
      flag[curr_index]= PROCESSED_BORDER_AIR;
      //cout<<"Processing_Border_Air_List - ending point \n";
}

void Navier::Ordered_New_Fluid_Cell_Initializer(int index, int curr_antecessors)
{ 
    // cout<<"Ordered_New_Fluid_Cell_Initializer - starting point \n";
    
    // cout<<"index ="<<index<<"\n";
     if(index ==  144856) cout<<"*****curr_antecessors = " << curr_antecessors <<"\n";
     static int dM1[3] ={1,2,0};
     static int symmetrical[6] ={1,0,3,2,5,4};
     static int orientation_sign[2] = {-1,1};
     int delabs, delmod;     
     double input_area;
     //double UVW_in;
     double UVW[3];
     double UVWout[3][2]={{0.0,0.0},{0.0,0.0},{0.0,0.0}}; 
     double coefsum [3][2]={{1.0E-15,1.0E-15},{1.0E-15,1.0E-15},{1.0E-15,1.0E-15}};
     //double coef[3]={0.0,0.0,0.0};
     double input_face_area[3]={0.0,0.0,0.0};
     double output_face_area[3]={0.0,0.0,0.0};
     int orientation_suc[6];
     int direction_suc[6];
     int input_area_suc[6];
     int sucessors_number=0;
     
     int output_area_index = 1;
     /*if(index== 18686)
     {
        for(int ifor=0; ifor<6;ifor++)
	{
	  for(int jfor=0; jfor<3;jfor++) 
	  {
	      //cout<<"griduvw ="<<index<<","<<ifor<<","<<jfor<<","<<index+ddebug[ifor]<<","<<gridUVW[jfor][index+ddebug[ifor]]<<"\n";
	  }
	}
     }
        for(int jfor=0; jfor<3;jfor++) 
	 {
	     cout<<"griduvw ="<<index<<","<<gridUVW[jfor][index]<<"\n";
	     if(isnan(gridUVW[jfor][index])){ cout<<"Initializer_start"; exit(0);} 
	 }*/
     for(int direction=0; direction<3;direction++)
     {
	 delabs = - d123[direction];
         delmod = 0;
	 int direction1 = dM1[direction]; 
	 int direction2 = dM1[direction1];     
	 for(int orientation=0; orientation<2;orientation++)
         {
	        //cout<<"Ordered_New_Fluid_Cell_Initializer - top inner inner for\n";
                if(flag[index+delabs] >= PROCESSED_BORDER_AIR || ((gridUVW[direction][index + delmod]* delabs) < 0) )
                  {
	             UVW[direction]= gridUVW[direction][index + delmod];
		    // UVW[direction]= gridCellUVW[direction][index];
		     //UVW_in= gridUVW[direction][index + delmod];
                     if(UVW[direction] * delabs < mepsilonxh)
                     {
		        int index_abs = index+delabs;
	                //cout<< index+delmod<<","<<index_mod+d123[direction1]<<","<<index_mod+d123[direction2]<<"\n";
			//cout<<"Ordered_New_Fluid_Cell_Initializer - Determining UVW\n";
			/*UVW[direction1] = fmax(gridUVW[direction1][index_abs],0.0);
			if(gridUVW[direction1][index_abs+d123[direction1]]< 0.0)
			  UVW[direction1]+=gridUVW[direction1][index_abs+d123[direction1]];
			UVW[direction2] = fmax(gridUVW[direction2][index_abs],0.0);
			if(gridUVW[direction2][index_abs+d123[direction2]]< 0.0)
			  UVW[direction2]+=gridUVW[direction2][index_abs+d123[direction2]]; */ 
			int den = 0;
 			UVW[direction1] = 0.0;
			if(gridUVW[direction1][index_abs] > 0) { UVW[direction1] += gridUVW[direction1][index_abs]; den = 1;}
			if(gridUVW[direction1][index_abs+d123[direction1]] < 0) { UVW[direction1] += gridUVW[direction1][index_abs+d123[direction1]]; den++;}
			if(den>0) UVW[direction1] /= den;
			UVW[direction2] = 0.0; den = 0;
			if(gridUVW[direction2][index_abs] > 0) { UVW[direction2] += gridUVW[direction2][index_abs]; den = 1;}
			if(gridUVW[direction2][index_abs+d123[direction2]] < 0) { UVW[direction2] += gridUVW[direction2][index_abs+d123[direction2]]; den++;}
			if(den>0) UVW[direction2] /= den;
			//UVW[direction1] = 0.5*(gridUVW[direction1][index_abs]+ gridUVW[direction1][index_abs+d123[direction1]]);
			//UVW[direction1] = gridCellUVW[direction1][index];
 			//UVW[direction2] = 0.5*(gridUVW[direction2][index_abs]+ gridUVW[direction2][index_abs+d123[direction2]]);
			//UVW[direction2] = gridCellUVW[direction2][index];
			if(index== 144856)
                        {
                           
	                    cout<<"index,direction,UVW ="<<index<<","<<direction<<","<<UVW[direction]<<","<<UVW[direction1]<<","<<UVW[direction2]<<"\n";
	                    cout<<index+delabs<<","<<flag[index+delabs]<<","<<flag[index+dj]<<"\n";
			    cout<<"UVW ="<<gridUVW[0][index]<<","<<gridUVW[1][index]<<","<<gridUVW[2][index]<<","<<gridUVW[0][index+1]<<","<<gridUVW[1][index+dj]<<
			    gridUVW[2][index+dk]<<"\n";
                        }
			//cout<<index+delabs<<","<<flag[index+delabs]<<"\n";
		   	if(gridPhi[index+delabs] > Phi_AIR_FLUID_limit)
			{
			     //cout<< output_area_index <<"\n";
			     input_area = output_area[symmetrical[output_area_index]][index+delabs];
			     //cout<< input_area <<"\n";
			}
			else
			{
			   input_area = 1.0;
			}
			
		     }//end-if(UVW[direction] * delabs < mepsilonxh)
		     else
		     {
			 delabs = -delabs;
	                 delmod = d123[direction];
			 output_area_index --;
		         continue;
		     }
		  }
		  else
		  {
		     if(flag[index+delabs] > AIR)
	             {
			    
			 orientation_suc[sucessors_number] = orientation;
			 direction_suc[sucessors_number] = direction;
			 input_area_suc[sucessors_number]= output_area_index;
			 
			 if(true)
		         {  
		            //cout<<"index ="<<index<<"\n";
			    //cout<<"sucessors_number ="<<sucessors_number<<"\n";
		            //cout<<"orientation_suc[itfor] ="<<orientation_suc[sucessors_number]<<"\n";
		            //cout<<"direction_suc[itfor] ="<<direction_suc[sucessors_number]<<"\n";
		            //cout<<" input_area_suc[sucessors_number] ="<< input_area_suc[sucessors_number]<<"\n";
		         }
			 sucessors_number ++;
		     }
		     delabs = -delabs;
	             delmod = d123[direction];
		     output_area_index --;
		     continue;
		  }    
           
	        //cout<<"Ordered_New_Fluid_Cell_Initializer - first if end\n";
	        
		
		double UVWinitial[3]={0.0,0.0,0.0};
                double  UVWmax = 0.0;
		int kmax;
		
                for(int kfor=0;kfor<3;kfor++)
		{
		    if(UVW[kfor] <0) UVWinitial[kfor] =1.0;
	            if(UVWmax < fabs(UVW[kfor]) )
		    {
              	       UVWmax = fabs(UVW[kfor]);
                       kmax  = kfor; 
	            }
                }
		//cout<<"Ordered_New_Fluid_Cell_Initializer - inner inner for 1 end\n";
		double other_side[3];
		 UVW[direction] = fabs(UVW[direction]);
		 other_side[direction] = UVW[direction]/UVWmax;
	       
                for(int kfor=direction1;kfor!=direction;kfor= (kfor+1)%3)
		{
                   other_side[kfor] = 1- fabs(UVW[kfor]/UVWmax);
                   
	        }
		if(kmax == direction)
                {
		   double input_full_face_area_0 = other_side[direction1]*other_side[direction2];
		   if(input_area <= input_full_face_area_0)
		   {
		   output_face_area[direction] = input_area ;
		   output_face_area[direction1] = 0.0;
		   output_face_area[direction2] = 0.0;
		   input_face_area[direction] = input_area;
		   input_face_area[direction1] = 0.0;
		   input_face_area[direction2] = 0.0;
		   }	
		   else 
		   {
		        input_face_area[direction] = input_full_face_area_0;
		        output_face_area[direction] = input_full_face_area_0;
		        double area_opposite =  (1 - other_side[direction1])*(1 - other_side[direction2]);
		        double a =(1 - input_full_face_area_0 + area_opposite)/(2.0 * area_opposite);
			//if(a > area_opposite) {a =area_opposite; cout<<" a > 1 \n"; int aa; cin>>aa;}
			if(a < 0.0) {a =0.0; cout<<" a < 0 \n"; int aa; cin>>aa;}
			double b= 1.0;
			if(fabs(area_opposite)>1E-16)
			 { 
			    if(input_area > 1.0) {input_area =1.0; cout<<" input_area > 1 \n"; int aa; cin>>aa;}
			    //if(input_area < 0.0) {input_area =0.0; cout<<" input_area < 0 \n"; /*int aa; cin>>aa;*/}
			    //if(area_opposite < 0.0) {cout<<" area_opposite < 0.0 < 0 \n"; int aa; cin>>aa;}
			    double delta = a*a - ((1.0 - input_area)/area_opposite);
			    if(delta >= 0) b = a - sqrt(delta);
			    else b = a;
			}
			else 
			{
			    if((1 - other_side[direction1]) > 1E-8) b = (1.0 - input_area)/(1 - other_side[direction1]);
			    else
			    {
			       if((1 - other_side[direction2]) > 1E-8) b = (1.0 - input_area)/(1 - other_side[direction2]);
			    }
			}    
			if(b > 1.0) {b =1.0; cout<<" b > 1 \n"; int aa; cin>>aa;} 
			if(b < 0.0) {b =0.0; cout<<" b < 0 \n"; int aa; cin>>aa;}      
		        //double b = a - sqrt(a*a - (1.0 - input_area)/area_opposite);
			//cout<<" I passed here: input_face_area[direction1] \n";
                        input_face_area[direction1] = 0.5 *(1-b)*((1 +  other_side[direction2])*(1-other_side[direction1]) - b *area_opposite);
			if( input_face_area[direction1] < 0.0) {input_face_area[direction1] =0.0; cout<<" input_face_area[direction1] < 0 \n"; int aa; cin>>aa;} 
			input_face_area[direction2] = fmax(0.0,input_area - input_full_face_area_0 - input_face_area[direction1]);
		        output_face_area[direction1] =  input_face_area[direction1]/ (1-other_side[direction1]);
		        output_face_area[direction2] = input_face_area[direction2]/(1-other_side[direction2]);
		   
		        //input_face_area[direction1] = 0.5 *(1 +  other_side[direction2])*(1-other_side[direction1]);
		        //input_face_area[direction2] = 0.5 *(1 +  other_side[direction1])*(1-other_side[direction2]);
		   }
		}
                else
		{	    
                   output_face_area[direction] = 0.0;
		   input_face_area[direction] = 0.0;
		   if(kmax == direction1)
		   {
		      double input_full_face_area_1 = 0.5 *(1 +  other_side[direction2]);
		      double output_full_face_area_1 = input_full_face_area_1 * other_side[direction];// = 0.5 *(1 +  other_side[direction2])*other_side[direction];
		      if(input_area <= input_full_face_area_1)
		      {
		          input_face_area[direction1]=  input_area;
		          input_face_area[direction2]=  0.0;  
		          output_face_area[direction1] = (input_area/input_full_face_area_1)* output_full_face_area_1;
		          output_face_area[direction2] = 0.0;
		      }
		      else
		      {
		          input_face_area[direction1]=  input_full_face_area_1;
		          input_face_area[direction2]=  input_area - input_full_face_area_1;
			  double output_full_face_area_2 = 0.5 * other_side[direction];
		          output_face_area[direction1] = output_full_face_area_1;
		          output_face_area[direction2] = (input_area - input_full_face_area_1)/(1.0 - input_full_face_area_1)* output_full_face_area_2;;
		      }
		   }
		   else
		   {
		      double input_full_face_area_2 = 0.5 *(1 +  other_side[direction1]);
		      double output_full_face_area_2 = input_full_face_area_2 * other_side[direction];// = 0.5 *(1 +  other_side[direction1])*other_side[direction];
		      if(input_area <= input_full_face_area_2)
		      {
		          input_face_area[direction2]=  input_area;
		          input_face_area[direction1]=  0.0;  
		          output_face_area[direction2] = (input_area/input_full_face_area_2)* output_full_face_area_2;
		          output_face_area[direction1] = 0.0;
		      }
		      else
		      {
		          input_face_area[direction2]=  input_full_face_area_2;
		          input_face_area[direction1]=  input_area - input_full_face_area_2;
			  double output_full_face_area_1 = 0.5 * other_side[direction];
		          output_face_area[direction2] = output_full_face_area_2;
		          output_face_area[direction1] = (input_area - input_full_face_area_2)/(1.0 - input_full_face_area_2)* output_full_face_area_1;;
		      }
		     
		     
		   }
		   /*{  
		      output_face_area[direction2] = 0.5 *(1 +  other_side[direction1])*other_side[direction];
		      output_face_area[direction1] = 0.5 * other_side[direction];
		      input_face_area[direction2] = 0.5 *(1 +  other_side[direction1]);
		      input_face_area[direction1] = 0.5 *(1 - other_side[direction1]);
		   }*/
		}
		//if(index == 18686)
                        //{
                           
	                    //cout<<"output face ="<<output_face_area[0]<<","<<output_face_area[1]<<","<<output_face_area[2]<<"\n";
	                    //cout<<"input face ="<<input_face_area[0]<<","<<input_face_area[1]<<","<<input_face_area[2]<<"\n";
			    //cout<<"index="<<index<<"\n";
                        //}
			
		//cout<<"Ordered_New_Fluid_Cell_Initializer - inner inner for 2 end\n";
		double input_area_sum =0.0;
		for(int lfor=0;lfor<3;lfor++) input_area_sum += input_face_area[lfor];
		if(input_area_sum > input_area)
		{
		   for(int lfor=0;lfor<3;lfor++)
                  {   
		      input_face_area[lfor]= input_face_area[lfor]* input_area/input_area_sum;
	              //cout<<"lfor-coef[lfor], UVWout,coefsum,="<<lfor<<","<<","<<UVWout[lfor][1-(int)UVWinitial[lfor]]<<","<<coefsum[lfor][1-(int)UVWinitial[lfor]]<<","<<"\n";
	              //cout<<"UVwmax,kmax,UVWinitial, otherside, direction="<<UVWmax<<","<<kmax<<","<<UVWinitial[lfor]<<","<<other_side[lfor]<<","<<direction<<"\n";
		  }
	          //cout<<" direction1,2,input_face_area="<<direction1<<","<<direction2<<","<<input_face_area[0]<<"\n";
	         // cout<<"gridPhi[index+delabs],flag[index+delabs] ="<<gridPhi[index+delabs]<<","<<flag[index+delabs]<<"\n";
		  //cout<<"input_area_sum,input_area ="<<input_area_sum<<","<<input_area<<"\n";
		  //int a; cin>>a;
                }
                for(int lfor=0;lfor<3;lfor++)
		{
		  int orientation_temp = 1 - (int)UVWinitial[lfor];
		  //coef[lfor] = input_area * input_face_area[lfor];
		  //UVWout[lfor][orientation_temp] += orientation_sign[orientation_temp]* UVW[direction]*coef[lfor];
		  if(input_face_area[lfor]< 0.0) cout<<"input_face_area negative="<< input_face_area[lfor]<<"\n";
		  UVWout[lfor][orientation_temp] += orientation_sign[orientation_temp]* UVW[direction]*input_face_area[lfor];
		  
		  
		  //UVWout[lfor][orientation_temp] += orientation_sign[orientation_temp]* UVW_in*input_face_area[lfor];
		  
		  coefsum[lfor][orientation_temp] += output_face_area[lfor];
		  
		  if(index== 144856)
                        {
                           
	                   cout<<" UVWout,coefsum,="<<lfor<<","<<orientation_temp<<","<<UVWout[lfor][orientation_temp]<<","<<coefsum[lfor][orientation_temp]<<"\n";
	                   cout<<"UVW[direction]*input_face_area[lfor]"<<UVW[direction]<<","<<input_face_area[lfor]<<"\n";
			    //cout<<" direction1,2,input_face_area="<<direction1<<","<<direction2<<","<<input_face_area[0]<<"\n";
			    //cout<<"gridPhi[index+delabs],flag[index+delabs] ="<<gridPhi[index+delabs]<<","<<flag[index+delabs]<<"\n";
                        }
		}
	        //if(input_face_area[0]<= 0.99)
		//{
		   //int x,y,z;
		  // GIJK(index,x,y,z);
		   //cout<<" index,x,y,z="<<x<<","<<y<<","<<z<<"\n";
		   //GIJK(index+delabs,x,y,z);
		   //cout<<" index+delabs,x,y,z="<<index+delabs<<","<<x<<","<<y<<","<<z<<"\n";
		   //cout<<" index,x,y,z="<<gridUVW[0][index+delabs]<<","<<gridUVW[1][index+delabs]<<","<<gridUVW[2][index+delabs]<<"\n";
		   //exit(1);
		 //}
	        delabs = -delabs;
	        delmod = d123[direction];
		output_area_index --;
		//cout<<"Ordered_New_Fluid_Cell_Initializer - inner inner for 3 end\n";
	  }// end for orientation
          output_area_index +=4;
     }//end for direction
     //cout<<"Ordered_New_Fluid_Cell_Initializer - Outer_for_end \n";
     for(int direction=0; direction<3;direction++)
     {
        delabs = - d123[direction];
        delmod= 0;
     	for(int orientation=0;orientation<2;orientation++)
     	{
                  if(flag[index+delabs] < PROCESSED_BORDER_AIR)
                  {
                     if(coefsum[direction][orientation]<0.0) cout<<"coefsum negative \n";		     
		     if(coefsum[direction][orientation]>epsilon)
		     {
		          if(gridUVW[direction] [index+delabs] < Phi_AIR_FLUID_limit || flag[index+delabs] == AIR)
			  {
	                     gridUVW[direction][index+delmod] = UVWout[direction][orientation];// /fmin(1.0,coefsum[direction][orientation]);
			     coefsum[direction][orientation] = 1.0;
			  }
			  else
			      gridUVW[direction][index+delmod] = UVWout[direction][orientation]/fmin(1.0,coefsum[direction][orientation]);
		     }
		  }
		     
		    /* if(direction == 0 && fabs(gridUVW[direction][index+delmod] - 0.3)> 0.001)
		     {
		       int x,y,z;
		       GIJK(x,y,z,index);
		       cout<< "Inicialization - griduvw  - 0="<<index<<","<<x<<","<<y<<","<<z<<","<<gridUVW[direction][index+delmod]<<"\n";
		       int a;
		       cin>>a;*/
		     //if(isnan(gridUVW[direction][index+delmod])){ cout<<"Initializer_end"; exit(0);} 
		   else
		   {
		      if(flag[index+delabs] >= FLUID && flag[index+delabs] < INFLOW )
		      {
		         if((gridUVW[direction][index+delmod] * UVWout[direction][orientation]) > 0)
			 {
			    if(fabs(gridUVW[direction][index+delmod]) > fabs(UVWout[direction][orientation]))
			    {
			       gridUVW[direction][index+delmod] = UVWout[direction][orientation]; 
			    } 
			 }
		       }   
                   }
		   if(index== 144856)
		   {
		       cout<<"d,o,UWWout"<<direction<<","<<orientation<<","<<UVWout[direction][orientation]<<"\n";
		       int x,y,z;
		       GIJK(index, x,y,z);
		       cout<< "Inicialization - griduvw  - 1="<<index+delmod<<","<<x<<","<<y<<","<<z<<","<<gridUVW[direction][index+delmod]<<","<<flag[index+delabs]<<"\n";
		       //int a;
		      // cin>>a;
		   }
                  
                  delabs = -delabs;
	          delmod = d123[direction];
     	}
	
	//cout<<"Ordered_New_Fluid_Cell_Initializer - inner inner for 4 end\n";
     }
    /* for(int direction=0; direction<3;direction++)
     {
        delabs = - d123[direction];
        delmod= 0;
     	for(int orientation=0;orientation<2;orientation++)
     	{
                  if(flag[index+delabs] < PROCESSED_BORDER_AIR)
                  {
                     if(coefsum[direction][orientation]<0.0) cout<<"coefsum negative \n";		     
                     if(coefsum[direction][orientation]>epsilon)
		     {
		          if(gridUVW[direction] [index+delabs] < Phi_AIR_FLUID_limit)
			  {
	                     gridUVW[direction][index+delmod] = UVWout[direction][orientation];// /fmin(1.0,coefsum[direction][orientation]);
			     coefsum[direction][orientation] = 1.0;
			  }
			  else
			      gridUVW[direction][index+delmod] = UVWout[direction][orientation]/fmin(1.0,coefsum[direction][orientation]);
		     }
		     else
		     {
		          //cout<<"coefsumzero index ="<<index<<"\n";
			  //gridUVW[direction][index+delmod] = 0.0;
		     }*/
		   /* if(direction == 1 && fabs(gridUVW[direction][index+delmod])> 0.001)
		     {
		       int x,y,z;
		       GIJK(x,y,z,index);
		       cout<< "Inicialization - griduvw  - 1="<<index<<","<<x<<","<<y<<","<<z<<","<<gridUVW[direction][index+delmod]<<"\n";
		       int a;
		       cin>>a;
		      }
		     */
		    /* if(direction == 0 && fabs(gridUVW[direction][index+delmod] - 0.3)> 0.001)
		     {
		       int x,y,z;
		       GIJK(x,y,z,index);
		       cout<< "Inicialization - griduvw  - 0="<<index<<","<<x<<","<<y<<","<<z<<","<<gridUVW[direction][index+delmod]<<"\n";
		       int a;
		       cin>>a;
		      }*/
		    /* if(isnan(gridUVW[direction][index+delmod])){ cout<<"Initializer_end"; exit(0);} 
                  }
                  
                  delabs = -delabs;
	          delmod = d123[direction];
     	}
	
	//cout<<"Ordered_New_Fluid_Cell_Initializer - inner inner for 4 end\n";
     }*/
     if(index== true)
		     {
		        //cout<<"gridPhi[index] ="<<gridPhi[index]<<"\n";
		        //cout<<"Phi_AIR_FLUID_limit ="<<Phi_AIR_FLUID_limit<<"\n";
		        //cout<<"index ="<<index<<"\n";
		     }
     if(gridPhi[index] > Phi_AIR_FLUID_limit)
     {
         
	 
	     for(int itfor=0;itfor<sucessors_number;itfor++) 
	     {
                
			 
		 //if(true)
		     //{  //cout<<" input_area_suc[itfor] ="<< input_area_suc[itfor]<<"\n";
		        //cout<<"index ="<<index<<"\n";
			//cout<<"sucessors_number ="<<sucessors_number<<"\n";
		        //cout<<"orientation_suc[itfor] ="<<orientation_suc[itfor]<<"\n";
		        //cout<<"direction_suc[itfor] ="<<direction_suc[itfor]<<"\n";
		        //cout<<" output_area[input_area_suc[itfor]][index] ="<< output_area[input_area_suc[itfor]][index]<<"\n";
			
		     //}
		 output_area[input_area_suc[itfor]][index] = fmin(1.0,coefsum[direction_suc[itfor]][orientation_suc[itfor]]);
             }
	 //cout<<"Ordered_New_Fluid_Cell_Initializer - inner inner for 5 end\n";	 
        
     }
     
     flag[index] = PROCESSED_BORDER_AIR;  
     //if(index ==31789) {cout<<" outflow_cell taken as air 1\n"; int a; cin>>a;}
     //cout<<"Ordered_New_Fluid_Cell_Initializer -  end\n"; 
}
	

void Navier::DoIt_otim()
{  
   cout<<"Do it starting-point \n";
   for( int ifor=0; ifor < dj; ifor++)
       cout<<"flag["<< ifor<<"]="<<flag[8*dk + ifor*dj +1]<<"\n";
   //if((num_iter_Navier == 100) || (num_iter_Navier == 101)){ cout<<*Border_Air_Non_Diagonal_size <<"\n"; int a; cin>> a;}
   int index,index1,init;    
   timestepmethod=0;
   List_Border_Air->clear();
   List_Border_Fluid->clear();
   //cout<<"Extremes of List ext-air->begin,air->end,fluid->begin,fluid->end"<<"\n";
   //cout<<List_Border_Face_Air->size()<<"\n";
   //cout	<<List_Border_Face_Fluid->size()<<"\n";
   
   //cout<<"size Border_Air="<<List_Border_Air->size()<<" \n";
   int itend = List_Border_Face_size_m1;
   for(int ifor=0; ifor<=itend;ifor++)
    {
       index = (*List_Border_Face_Air)[ifor];
       index1 = (*List_Border_Face_Fluid)[ifor];
       init =0;
       if(flag[index] > AIR_FLUID_LIMIT)
       {
           cout<<"Air_labelled_Fluid = "<< index <<","<<flag[index]<<"\n";
	   cout<<"Air_labelled_Fluid = "<< index1 <<","<<flag[index1]<<"\n";
	   int x,y,z;
	   GIJK(index,x,y,z);
	   cout<<"x,y,z(index)="<<x<<","<<y<<","<<z<<"\n";
	   GIJK(index1,x,y,z);
	   cout<<"x,y,z(index1)="<<x<<","<<y<<","<<z<<"\n";
           int a; cin>>a;
       }
       if(flag[index]%incrementx2 < increment) 
       {
           flag[index] += increment;
           List_Border_Air->push_back(index);
           int delta_index = abs(index - index1);
           init=(delta_index==dk)?4:((delta_index==dj)?2:1);
           List_Cost_Air.push_back(init);
       }
        //cout<<"flag[index1]="<<flag[index1]<<" \n";
	//cout<<"flag[index1]% incrementx2 ="<<flag[index1]% incrementx2 <<" \n";
      if(flag[index1] % incrementx2 < increment) 
       {
           flag[index1] += increment;
           List_Border_Fluid->push_back(index1);
	   if(init == 0)
	   {
              int delta_index = abs(index - index1);
              init=(delta_index==dk)?4:((delta_index==dj)?2:1);
           }
	   List_Cost_Fluid.push_back(init);
       }
    }
    //cout<<"Do it First loop end \n";
    //cout<<"size Border_Air="<<List_Border_Air->size()<<" \n";
    //num_direct_neighbors_air = *List_Border_Air->size();
    *Border_Air_Non_Diagonal_size = (int) List_Border_Air->size();
     //cout<<"size Border_Fluid ="<<List_Border_Fluid->size()<<" \n";
    *Border_Fluid_Non_Diagonal_size = (int) List_Border_Fluid->size();  	 
    Diagonal_Border_Cells(List_Border_Air,List_Cost_Air,List_AirToFluid,1);
    //num_direct_neighbors_fluid = List_Border_Fluid->size();
    Diagonal_Border_Cells(List_Border_Fluid,List_Cost_Fluid,List_FluidToAir,-1);
    //cout<<"size Cost_Air="<<List_Cost_Air.size()<<" \n";
    //cout<<"size Border_Air="<<List_Border_Air->size()<<" \n";
    //cout<<"size Cost_Fluid="<<List_Cost_Fluid.size()<<" \n";
    //cout<<"size Border_Fluid="<<List_Border_Fluid->size()<<" \n";
   //cout<<"Diagonal Border Cells called \n";
   for(int itfor=0; itfor<= List_Border_Face_size_m1; itfor++)
   {
       bool flag_delete = false;
       int index0 = (*List_Border_Face_Air)[itfor];
       int index1 = (*List_Border_Face_Fluid)[itfor];
       if((gridPhi[index0]<0))
       {
           if(flag[index0]< AIR_FLUID_LIMIT)
	   {
	       if((flag[index0]%incrementx4) < incrementx2)
	       {
	           List_AirToFluid.push_back(index0);
	           flag[index0] += incrementx2;
	       }
	       flag_delete = true;
	    }
       }
       //cout<<"List AirToFluid Air_Part \n";
       if(gridPhi[index1]>= 0)
       {
	   if(flag[index1]> AIR_FLUID_LIMIT)
	   {
	       if(flag[index1] % incrementx4 < incrementx2)
	       {
	       List_FluidToAir.push_back(index1);
	       flag[index1] += incrementx2;
	       }
	       if(flag_delete)
	       {
	           int index_temp = (*List_Border_Face_Air)[itfor];
		       (*List_Border_Face_Air)[itfor]=(*List_Border_Face_Fluid)[itfor];
		       (*List_Border_Face_Fluid)[itfor]=index_temp;
		       flag_delete = false;
	       }
	       else
	           flag_delete = true;
	   }
       }
       if(flag_delete)
       {
         (*List_Border_Face_Air)[itfor]= (*List_Border_Face_Air)[List_Border_Face_size_m1];
         (*List_Border_Face_Fluid)[itfor]= (*List_Border_Face_Fluid)[List_Border_Face_size_m1];
         List_Border_Face_size_m1--;
	 //cout<<" List_Border_Face_size_m1 ="<< List_Border_Face_size_m1<<"\n";
	 itfor--; 
       }
       
    } 
   //cout<<"List Fluid_to_Air processed \n";
   itend = List_FluidToAir.size();
   for(int itfor=0;itfor< itend;itfor++)
   {
       int index = List_FluidToAir[itfor];
       {
	   for(int jt=0;jt<=2;jt++)
           {
              int index1 = index + d123[jt];
              for(int kt=1;kt<=2;kt++)
              {    
                 if((gridPhi[index1]<= 0) && (flag[index1]>= AIR)) 
                 {
		    List_Border_Face_size_m1++;
                    if(List_Border_Face_size_m1 >= (int)List_Border_Face_Air->size())
		    {
		        List_Border_Face_Air->push_back(index);
                        List_Border_Face_Fluid->push_back(index1);
		    }
		    else
		    {
		        (*List_Border_Face_Air)[List_Border_Face_size_m1] = index;
                        (*List_Border_Face_Fluid)[List_Border_Face_size_m1]= index1;
		    }
		      
                 }
                 index1 = index - d123[jt];
              }    
           }
        }
    }  
   itend = List_AirToFluid.size();
   for(int itfor=0;itfor<itend;itfor++)
   {
       int index = List_AirToFluid[itfor];
       {
	   for(int jt=0;jt<=2;jt++)
           {
              int index1 = index + d123[jt];
              for(int kt=1;kt<=2;kt++)
              {    
                 if((gridPhi[index1]> 0) && (flag[index1]>= AIR)) 
                 {
		    /* if(flag[index1] > AIR_FLUID_LIMIT)
                    {
                        cout<<"Air_labelled_Fluid = "<< index <<","<<flag[index]<<"\n";
	                cout<<"Air_labelled_Fluid = "<< index1 <<","<<flag[index1]<<"\n";
	                int x,y,z;
	                GIJK(index,x,y,z);
	                cout<<"x,y,z(index)="<<x<<","<<y<<","<<z<<"\n";
	                GIJK(index1,x,y,z);
	                cout<<"x,y,z(index1)="<<x<<","<<y<<","<<z<<"\n";
                        int a; cin>>a;
                    }*/
                    List_Border_Face_size_m1++;
                    if(List_Border_Face_size_m1 >= (int)List_Border_Face_Air->size())
		    {
		        List_Border_Face_Air->push_back(index1);
                        List_Border_Face_Fluid->push_back(index);
		    }
		    else
		    {
		        (*List_Border_Face_Air)[List_Border_Face_size_m1] = index1;
                        (*List_Border_Face_Fluid)[List_Border_Face_size_m1]= index;
		    }
                 }
                 index1 = index - d123[jt];
              }    
           }
        }
    } 
   double UVWmax_newFluidlist= 0,gridmax[6];
   int index_max = dj+dk+1;
   for(int itfor=0; itfor<List_Fluid_end_zero;itfor++)
   {
       index=(*List_Fluid_zero)[itfor];
       for(int jt=0;jt<=2;jt++)
       {
           if(fabs(gridUVW[jt][index])>UVWmax_newFluidlist)
	   {
	         UVWmax_newFluidlist= fabs(gridUVW[jt][index]);
	         index_max =index;
	         gridmax[0]= gridUVW[0][index];
	         gridmax[1]= gridUVW[1][index];
	         gridmax[2]= gridUVW[2][index];
	         gridmax[3]= gridUVW[0][index+1];
	         gridmax[4]= gridUVW[1][index+dj];
	         gridmax[5]= gridUVW[2][index+dk];
	   }
       }
   }
   cout<<"UVWmax_newFluidlist ="<<UVWmax_newFluidlist<<"\n";  
   for(int ifor = 0; ifor < 6; ifor++)  cout<<"gridmax ="<<gridmax[ifor]<<"\n";
   int x,y,z;
   GIJK(index_max,x,y,z);
   cout<<"x,y,z="<<index_max<<","<<x<<","<<y<<","<<z<<"\n";
   cout<<"flags="<<flag[index_max-1]<<","<<flag[index_max+1]<<","<<flag[index_max-dj]<<","<<flag[index_max+dj]<<","<<flag[index_max-dk]<<","<<flag[index_max+dk]<<"\n";
    
   Processing_Border_Air_List();     
   //cout<<"\n Doit Sai Processing_Border_Air \n";   
   //if(num_iter_Navier == 2) exit(0);        
  
   //cout<<"\n Doit Sai Call_CompFGH \n";  
   
   //cout<<"List_Border_Air->size()="<<List_Border_Air->size()<<"\n";
   //cout<<"List_Border_Fluid->size()="<<List_Border_Fluid->size()<<"\n";
   itend =(int)List_Border_Air->size();
   //cout<<" List_Border_Air \n";
   for(int itfor=0; itfor< itend; itfor++)
   {
      index = (*List_Border_Air)[itfor]; 
      //cout<<"index_Air="<<itfor<<","<<index<<"\n";
      /*int x,y,z;
        GIJK(index,x,y,z);
        //if(y>=14 && y<=18 && z>=14 && z<=18)
	if(true)
       {
          for(int jt=0;jt<=2;jt++)
          {
            
	            cout<<"x,y,z,grid ="<<x<<","<<y<<","<<z<<","<<gridUVW[jt][index]<<"\n";
	            if(isnan(gridUVW[jt][index])){ cout<<"Processing_Border_Air"; exit(0);}              
            } 
	    
        } */
      //sleep(1000);
     
      if(flag[index] <  SLIP_AIR)
      {
          //if(index ==31789) {cout<<" outflow_cell taken as air 2\n"; int a; cin>>a;}
	 flag[index] = AIR;
      }
      else
      {
          int mask = flag[index]% increment;
          flag[index] -= ((flag[index]% incrementx4) -mask);
       	 //flag[index] = flag[index]& erase_mask; 
      }
   }

   //cout<<" List_Border_Fluid \n";
   itend =(int)List_Border_Fluid->size();
   for(int itfor=0; itfor< itend; itfor++)
   {
       index = (*List_Border_Fluid)[itfor];
      /*  int x,y,z;
        GIJK(index,x,y,z);
       // if(y>=14 && y<=18 && z>=14 && z<=18) 
       if(true)
       {
          for(int jt=0;jt<=2;jt++)
          {
            
	            cout<<"x,y,z,grid ="<<x<<","<<y<<","<<z<<","<<gridUVW[jt][index]<<"\n";
	            if(isnan(gridUVW[jt][index])){ cout<<"Processing_Border_Air"; exit(0);}              
            } 
	    
        } */
       //cout<<"index_Fluid="<<itfor<<","<<index<<"\n";
       int mask = flag[index]% increment;
       flag[index] -= ((flag[index]% incrementx4) -mask);
       //flag[index] = flag[index]&erase_mask;; 
   }
    //cout<< "\n Doit Flags Corrected\n";
   //cout<<"List_FluidToAir.size()="<<List_FluidToAir.size()<<"\n"; 
   //cout<<"List_AirToFluid.size()="<<List_AirToFluid.size()<<"\n";
  // CompRHS();// Adams-Bashfort or Runge-Kutta Compute right hand side for poisson equation
  //Call_CompFGH();
   for(unsigned int itfor=0; itfor<List_FluidToAir.size(); itfor++) 
   {
      index = List_FluidToAir[itfor]; 
      if( flag[index] > AIR_FLUID_LIMIT)
      {   
          int maskx2 = 2*(flag[index]% increment);   
          flag[index] = 0x10000 - flag[index]+ maskx2;
          P[index] = 0.0;
      }
   }
   itend = List_AirToFluid.size();
   int List_Fluid_size =(int)List_Fluid_zero->size();
   //cout<<"List_Fluid_end_zero="<<List_Fluid_end_zero<<"\n";
   int itnewfluid=0;
   for(int itfor=0; itfor< itend; itfor++)
   {
      index = List_AirToFluid[itfor];
       if(flag[index] < AIR_FLUID_LIMIT)
       { 
          int maskx2 = 2*(flag[index]% increment);
          flag[index] = 0x10000 - flag[index] + maskx2;
	  if(flag[index] < FLUID)
	   {
	      List_Slip_Outflow_Fluid.push_back(index);
	     
	     if(flag[index]<SLIP_FLUID)
	     { 
	        //cout<<"List_Outflow_Fluid.push_back_antes\n";
	        List_Outflow_Fluid.push_back(index);
	       //cout<<"List_Outflow_Fluid.push_back_depois\n";
	     }
	   }
          int itfluid = List_Fluid_end_zero + itnewfluid;
	  itnewfluid ++;
          if(itfluid >=List_Fluid_size)
             List_Fluid_zero->push_back(index);
          else (*List_Fluid_zero)[itfluid] = index;
	}
   }
   cout<<"num_iter_Navier ="<<num_iter_Navier<<"\n";
   List_Fluid_end_zero += itnewfluid;
   //cout<<"List_Fluid_end_zero="<<List_Fluid_end_zero<<"\n";
   
   UVWmax_newFluidlist = 0.0;
   index_max=dj+dk+1;
   for(int itfor=0; itfor<List_Fluid_end_zero;itfor++)
   {
       index=(*List_Fluid_zero)[itfor];
       for(int jt=0;jt<=2;jt++)
       {
           if(fabs(gridUVW[jt][index])>UVWmax_newFluidlist)
	   {
	         UVWmax_newFluidlist= fabs(gridUVW[jt][index]);
	         index_max =index;
	         gridmax[0]= gridUVW[0][index];
	         gridmax[1]= gridUVW[1][index];
	         gridmax[2]= gridUVW[2][index];
	         gridmax[3]= gridUVW[0][index+1];
	         gridmax[4]= gridUVW[1][index+dj];
	         gridmax[5]= gridUVW[2][index+dk];
	   }
       }
   }
   cout<<"UVWmax_newFluidlist ="<<UVWmax_newFluidlist<<"\n";  
   for(int ifor = 0; ifor < 6; ifor++)  cout<<"gridmax ="<<gridmax[ifor]<<"\n";
   GIJK(index_max,x,y,z);
   cout<<"x,y,z="<<index_max<<","<<x<<","<<y<<","<<z<<"\n";
   cout<<"flags="<<flag[index_max-1]<<","<<flag[index_max+1]<<","<<flag[index_max-dj]<<","<<flag[index_max+dj]<<","<<flag[index_max-dk]<<","<<flag[index_max+dk]<<"\n"; 
    //cout<<"num_iter_Navier ="<<num_iter_Navier<<"\n";  
    //cout<< "\n Doit Atualizei List_Fluid_zero \n";
  //CompRHS();// Adams-Bashfort or Runge-Kutta Compute right hand side for poisson equation
  Call_CompFGH();//calculate temporary velocity fields (in space) due to chorin projection method (see User's Guide)
  
  //cout<<"\n Doit -  Sai Call_CompFGH() \n"; 
  /* for(int itfor=0; itfor<List_Fluid_end_zero;itfor++)
   {
       index=(*List_Fluid_zero)[itfor];
       for(int jt=0;jt<=2;jt++)
       {
           cout<<itfor<<","<<index<<","<<gridUVW[jt][index]<<","<<flag[index]<<"\n";
	   if(abs(gridUVW[jt][index]) > 10.0) exit(1);
       }
   } */
   Call_Update_Slip_Outflow_Cells();
   //cout<<"\n Doit -  Sai Call_Update_Slip_Outflow_Cells \n";  
   update_pressure = true;
   CompRHS();// Adams-Bashfort or Runge-Kutta Compute right hand side for poisson equation
   
   //if(num_iter_Navier>2000) exit(1);
  // cout<<"\n Doit -  Sai RHS \n";	 
   Poisson(); // Solve poisson equation 
   //cout<<"\n Doit - Sai Poisson \n"; 	   
   AdapUVW_Time_step(); // Compute new velocity (velocity correction by pressure-gradient) 
   //cout<<"\n Doit - Sai AdapUVW \n";
   update_pressure = false;
   CompRHS();
  for(unsigned int itfor=0; itfor< List_Outflow_Fluid.size(); itfor++) 
   {
      cout<<" List_Outflow_Fluid.size()"<< List_Outflow_Fluid.size()<< "\n";
      index =  List_Outflow_Fluid[itfor];
      if(flag[index]%incrementby2<incrementby4)
      {
         int *initial_pointer_table,*pointer_table;
         initial_pointer_table = &(neighbors_table[flag[index]%incrementby2][0]);
         for(pointer_table = initial_pointer_table + 1; pointer_table<=initial_pointer_table + *initial_pointer_table;pointer_table++)
         {
            int neighbor_temp = *pointer_table;
	    int index_neighbor_temp = index + cell_contact_face_delta[neighbor_temp];
	    if( flag[index_neighbor_temp] >= FLUID)
	      *Total_outflow += fabs(gridUVW[neighbor_temp/2][index_neighbor_temp])*dt;
	 }
       }
       else
       {
            int neighbor_temp = flag[index]%incrementby4;
            int index_neighbor_temp = index + cell_contact_face_delta[neighbor_temp];
	    if( flag[index_neighbor_temp] >= FLUID)
	      *Total_outflow += fabs(gridUVW[neighbor_temp/2][index_neighbor_temp])*dt;
       }
    }      
	      
	  
  S->delt_old=S->delt; // old time step needed for adams-bashfort          
  timestepmethod|=1;//used to mark the first time step for AB discretization and data output 
  S->timesteps++; //tag which defines if F,G,H arrays are from actual or previous time step  
  *List_Fluid_end = List_Fluid_end_zero;
  List_AirToFluid.clear(); // clear= erase(List.begin(),List.end()); 
  List_FluidToAir.clear();
  List_Cost_Air.clear();
  List_Cost_Fluid.clear();
  List_Border_Slip_Outflow_Air.clear(); 
  if(gridPhi == gridPhi1) gridPhi = gridPhi2;
  else gridPhi = gridPhi1;
  num_iter_Navier++;
  //cout<<"Do it end ="<<List_Fluid_end_zero<<"\n";
}
//
void Navier::Diagonal_Border_Cells(std::vector<int> *List,std::vector<int>& List_Cost,std::vector<int>& Target_List, int sign)
{   
     //cout<<"Diagonal_Border_Cells -starting point \n";
    
    static int bit[3]= {1,2,4};
    unsigned int inifor = 0;
    int endfor = (int)List->size();
   
    for(int ht=0; ht<=1; ht++)
    {
       for(int itfor=inifor;itfor<endfor;itfor++)
       {   
           int index = (*List)[itfor];
           int curr_seq_bit = List_Cost[itfor];
           for(int jt=2;jt>=0;jt--)
           {
              if ((curr_seq_bit/bit[jt]) == 0) 
              {
                  int index1 = index+d123[jt];
                  for(int kt=1;kt<=2;kt++)
                  {   
                     if(flag[index1]%incrementx2 < increment)
		     {
		         if((AIR_FLUID_LIMIT - flag[index1])*sign > 0) 
                         {
                            if((gridPhi[index1])*sign < 0)
			    {
			       if(flag[index1] % incrementx4 < incrementx2)
			       {
			           Target_List.push_back(index1);
			           flag[index1] += incrementx2;
			       }
			    }
			    List->push_back(index1);
                            List_Cost.push_back(List_Cost[itfor] + bit[jt]); 
                            flag[index1] += increment;  
                         }
                     }
		     index1 = index - d123[jt];
                  }  
	       }
	       curr_seq_bit = curr_seq_bit % bit[jt];
           }
       }
       inifor = endfor;
       endfor = (int)List->size();
   }
  // cout<<"Diagonal_Border_Cells -ending point \n";
}
//Compute auxiliary-velocities without free-boundaries ##
void Navier::Call_CompFGH()//calculate temporary velocity fields (in space) due to chorin projection method (see User's Guide) 
{ 
   cout<<" Call_CompFGH - starting point"<<"\n";
   //cout<<"List_Fluid_end_zero ="<<List_Fluid_end_zero<<"\n";
       
   Der_Conv_Max=0.0;
   for(int itfor=0; itfor<List_Fluid_end_zero;itfor++)
   { 
       int index=(*List_Fluid_zero)[itfor];
       //cout<<"index ="<<index<<"\n";
       //cout<<"flag[index] ="<<flag[index]<<"\n";
       if((flag[index]>=FLUID) && (flag[index]<INFLOW))             
       {
		CompFGH(index); 
		for(int i0=0; i0<=2; i0++)
		{
		    if(flag[index + d123[i0]]< SLIP_FLUID)
		    {
		       init_FGH=end_FGH=i0;
		       CompFGH(index + d123[i0]);
		     }
		 }
		 init_FGH=0; end_FGH=2;
	         
	}
	else
	{
	   
	   List_Fluid_end_zero = List_Fluid_end_zero -1;
	   //cout<<"List_Fluid_end_zero COMPFGH="<<List_Fluid_end_zero<<"\n";
	   (*List_Fluid_zero)[itfor]=(*List_Fluid_zero)[List_Fluid_end_zero];
	   //cout<<"(*List_Fluid)[itfor]="<< (*List_Fluid)[itfor]<<"\n";
	   itfor--;
	} 
	 
    }
    cout<<" Call_CompFGH - ending point"<<"\n";
    for(int kt=0;kt<=2;kt++)
    {
      gridUVW[kt] += deltaUVW;
    }
    deltaUVW = -deltaUVW;
    cout<<" Der_Conv_Max= "<< Der_Conv_Max<<" \n";
    bool flag_temp = true;
    for(int ifor =0; ifor<=6; ifor++)
    {
       if( fabs(Der_Conv_Max) <= fabs(grid_max_advected[ifor])) flag_temp = false;
	cout<<"grid_max_advected["<<ifor<<"] ="<< grid_max_advected[ifor]<<"\n";
	cout<<"flag_max_advected["<<ifor<<"] ="<< flag_max_advected[ifor]<<"\n";
    }
    if(flag_temp){ int a; cin>>a;}	     
}
void Navier::CompFGH(int index)
{  
	//cout<<" CompFGH - starting point"<<"\n";
	//static int hx2 = 2*h;
        S->nu = 0.0;
	static int dx2123[3] ={2,2*dj,2*dk};
	static int iM1[3] ={1,2,0};
	
	
	flag_border_advection =true;  
	for(int i0=init_FGH; i0<=end_FGH; i0++)
	{  
	   Sum_coef = 0.0;
           int d0 =d123[i0];
	   if(flag[index-d0] == INFLOW) continue;
	   int d20=2*d0;
	           Der_Conv = 0;
		   L = gridUVW[i0][index-d0];
		   M = gridUVW[i0][index];
		   R = gridUVW[i0][index+d0];
		   if((flag[index-d20]>= FLUID) && (flag[index+d0] >= FLUID))
		   {
		      LL=gridUVW[i0][index-d20] ;
		      RR=gridUVW[i0][index+d20] ;
		      //flag_border_advection = false;
		   }
		   else 
		   {  
		      //if(flag[index-d0] == INFLOW) continue;
		      flag_border_advection = true;  
		   }
		   Laplacian = (L+R-2*M)/h2; 
		     KLL = L;
		     KRR = R;
		     KM = M;
		     /*KL=(L+M)/2;
		     KR=(R+M)/2;*/
		     L=KL=(L+M)/2;
		     R=KR=(R+M)/2;
		   //cout<<"cout<<L<<<<R<<<<M<<<<KL<<<<KR<<<<Laplacian<<<<index<<\n";
		   //cout<<L<<","<<R<<","<<M<<","<<KL<<","<<KR<<","<<Laplacian<<","<<index<<"\n";
		   
		   Comp_Der_Conv();
		   //
		   int I12 = iM1[i0];
		   do
		   {
			int D2122 = dx2123[I12];
			int D12 = d123[I12];
			/*L=gridUVW[i0][index-D12];
			R=gridUVW[i0][index+D12];*/
			L = 0.5 * (M + gridUVW[i0][index-D12]);
			R = 0.5 * (M + gridUVW[i0][index + D12]);
			if((flag[index-D2122] >= FLUID) && (flag[index+D2122] >= FLUID))
			{
				LL=gridUVW[i0][index-D2122];
				RR=gridUVW[i0][index+D2122];
				//flag_border_advection = false;
			 }
			 else  flag_border_advection = true;
			 
		   	Laplacian +=(L+R-2*M)/h2; 
			
		       /* KLL = 0.5*(gridUVW[I12][index - D12] + gridUVW[I12][index - d0 - D12]);
			KRR = 0.5 *(gridUVW[I12][index+ D12] + gridUVW[I12][index - d0 + D12]); 
		        KM = 0.5*(gridUVW[I12][index]+ gridUVW[I12][index-d0]);
			KL = 0.5*(KM + KLL);
			KR = 0.5*(KM + KRR);*/
			KL = 0.5*(gridUVW[I12][index]+ gridUVW[I12][index-d0]);
			KR = 0.5 *(gridUVW[I12][index+ D12] + gridUVW[I12][index - d0 + D12]);
			KM = 0.5 *(KL +KR);
			
			//cout<<"cout<<L<<<<R<<<<M<<<<KL<<<<KR<<<<Laplacian<<<<index<<\n";
		        //cout<<L<<","<<R<<","<<M<<","<<KL<<","<<KR<<","<<Laplacian<<","<<index<<"\n"; 
			 
			Comp_Der_Conv(); 
			I12 =iM1[I12];  
		    }
		    while(I12 != i0);
		  double gridadvected = (gridUVW[i0][index] - (Der_Conv * dt))/ (1 + ((Sum_coef/h)*dt));
		  
		  if(fabs(gridadvected)> fabs(Der_Conv_Max))
		  {
		     Der_Conv_Max = gridadvected;
		     grid_max_advected[0] = gridUVW[i0][index];
		     grid_max_advected[1] = gridUVW[i0][index - 1];
		     grid_max_advected[2] = gridUVW[i0][index + 1];
		     grid_max_advected[3] = gridUVW[i0][index - dj];
		     grid_max_advected[4] = gridUVW[i0][index + dj];
		     grid_max_advected[5] = gridUVW[i0][index - dk];
		     grid_max_advected[6] = gridUVW[i0][index + dk];
		     flag_max_advected[0] = flag[index];
		     flag_max_advected[1] = flag[index - 1];
		     flag_max_advected[2] = flag[index + 1];
		     flag_max_advected[3] = flag[index - dj];
		     flag_max_advected[4] = flag[index + dj];
		     flag_max_advected[5] = flag[index - dk];
		     flag_max_advected[6] = flag[index + dk];
		   }
		   
		   gridUVW[i0][index+deltaUVW] = (viscosity * Laplacian + Ext_Force[i0])* dt + gridadvected;
		  //gridUVW[i0][index+deltaUVW] = (viscosity * Laplacian - Der_Conv + Ext_Force[i0])* dt + gridUVW[i0][index];
		  
		  /*if(i0 == 1 && fabs(gridUVW[i0][index+deltaUVW])> 0.001)
		     {
		       cout<< "CompFGH - griduvw ="<<index<<","<<gridUVW[i0][index+deltaUVW]<<"\n";
		       int a;
		       cin>>a;
		      } 
		   if(i0 == 0 && fabs(gridUVW[i0][index+deltaUVW] -0.3)> 0.001)
		     {
		       cout<< "CompFGH - griduvw - 0 ="<<index<<","<<gridUVW[i0][index+deltaUVW]<<"\n";
		       int a;
		       cin>>a;
		      } */
		 // cout<<" gridUVW[i0][index+deltaUVW] ="<<gridUVW[i0][index+deltaUVW]<<"\n";                    
        }
	//cout<<" CompFGH - ending point"<<"\n";
	//cout<<"h,h2 ="<<h<<","<<h2<<"\n";
}
//
//
void  Navier::Comp_Der_Conv()
{
//####################################################################
//#                                                                  #
//# Discretization of convective Terms:                              #
//#                                                                  #
//#       e.g. discretization of \partial(uu)/\partial x             #
//#                                                                  #
//#    ------------- -------------- -------------- --------------    #
//#   |             |              |              |              |   #
//#   |             |              |              |              |   #
//#   |             |              |              |              |   #
//#   LL            L       KL     M      KR      R              RR  #
//#   |             |              |              |              |   #
//#   |             |              |              |              |   #
//#   |             |              |              |              |   #
//#    ------------- -------------- -------------- --------------    #
//#    \____  _____/ \_____  _____/ \______  ____/ \_____  _____/    #
//#         \/             \/              \/            \/          #
//#         DD             D               DP            PP          #
//#                                                                  #
//####################################################################
    

 // Sharp and Monotonic Algorithm for Realistic Transport (SMART) 2nd-Order
          
     if(flag_border_advection == false)
     {  
         if(KR>=0)
	 {       
            Phi_R=(M-L)/(((R-L)!=0.0)?(R-L):1E-15);
	    if(Phi_R>=0)
	    {
	       if(Phi_R<(3./74))
	       {
                  UR=10*M-9*L;
	       }
               else 
	       {
                   if(Phi_R<(5./6))
		   {
                      UR=(3.*R+6.*M-L)/8.; 
                     //UR=(1-UR)*((3./8)*R+(6./8)*M-(1./8)*L)+UR*((3./8)*M+(6./8)*R-(1./8)*RR);
		   }
                   else 
		   {
                      if(Phi_R<=1.0)
		      {
                         UR=R;
		      }
                      else 
                      {
		         UR=M;
		      }
		    }
		 }
	     }
	     else
	     {
	        UR=M;
	     }
	 }
	 else
	 {
	    Phi_R=(R-RR)/(((M-RR)!=0.0)?(M-RR):1E-15);
	    if(Phi_R>=0)
	    {
	       if(Phi_R<(3./74))
	       {
                  UR=10*R-9*RR;
	       }
               else 
	       {
                   if(Phi_R<(5./6))
		   {
                       UR= (3.*M+6.*R-RR)/8.;
                      //UR=(1-UR)*((3./8)*R+(6./8)*M-(1./8)*L)+UR*((3./8)*M+(6./8)*R-(1./8)*RR);
		   }
                   else 
		   {
                      if(Phi_R<=1.0)
		      {
                         UR=M;
		      }
                      else 
		      {
                         UR=R;
		      }
		   }
		}
	     }
	     else
	     {
	        UR=R;
	     }
	  }
	    
	 if(KL>=0)
	 {
	   
            Phi_L=(L-LL)/(((M-LL)!=0.0)?(M-LL):1E-15);
	    if(Phi_L>=0)
	    {
	        if(Phi_L<(3./74))
		{
                   UL=10*L-9*LL;
		}
                else
		{ 
                   if(Phi_L<(5./6))
		   {
                      UL=(3.*M+6.*L-LL)/8.;
                    //UL=(1-UL)*((3./8)*M+(6./8)*L-(1./8)*LL)+UL*((3./8)*L+(6./8)*M-(1./8)*R);
		   }
                   else 
		   {      
                      if(Phi_L<=1.0)
		      {
                         UL=M;
		      }
                      else 
		      {
                         UL=L;
		      }
		    }
		 }
	      }
	      else
	      {
	        UL=L;
	      }
	 }  
	 else
	 {
	    Phi_L=(M-R)/(((L-R)!=0.0)?(L-R):1E-15);
            if(Phi_L>=0)
	    {
	       if(Phi_L<(3./74))
	       {
                  UL=10*M-9*R;
	       }      
               else
	       { 
                  if(Phi_L<(5./6))
		  {
                     UL=(3.*L+6.*M-R)/8.;
                    //UL=(1-UL)*((3./8)*M+(6./8)*L-(1./8)*LL)+UL*((3./8)*L+(6./8)*M-(1./8)*R);
		  }
                  else 
		  {
                     if(Phi_L<=1.0)
		     {
                        UL=L;
		     }
                     else 
		     {
                        UL=M;
		     }
		  }
	       }
	    }
	    else
	    {
	      UL=M;
	    }
	  } 
	       
	
         Der_Conv += (KR*UR-KL*UL)/h;
	 Sum_coef = KL - KR;
     }
     /* Otiginal Version
     // Sharp and Monotonic Algorithm for Realistic Transport (SMART) 2nd-Order
     UR=(KR>=0) ? 0:1;
     UL=(KL>=0) ? 0:1;      
     if(flag_border_advection == false)
     {         
         Phi_R=(1-UR)*(M-L)/(((R-L)!=0.0)?(R-L):1E-15)  + UR*(R-RR)/(((M-RR)!=0.0)?(M-RR):1E-15);
         Phi_L=(1-UL)*(L-LL)/(((M-LL)!=0.0)?(M-LL):1E-15) + UL*(M-R)/(((L-R)!=0.0)?(L-R):1E-15);
      
         if((Phi_R>=0) && (Phi_R<(3./74)))
            UR=(1-UR)*(10*M-9*L)+UR*(10*R-9*RR);
         else 
            if((Phi_R>=(3./74)) && (Phi_R<(5./6)))
               UR=(1-UR)*((3.*R+6.*M-L)/8.) + UR*((3.*M+6.*R-RR)/8.);
               //UR=(1-UR)*((3./8)*R+(6./8)*M-(1./8)*L)+UR*((3./8)*M+(6./8)*R-(1./8)*RR);
            else 
               if((Phi_R>=(5./6)) && (Phi_R<=1))
                  UR=(1-UR)*R+UR*M;
               else 
                  UR=(1-UR)*M+UR*R;
	
          if((Phi_L>=0) && (Phi_L<(3./74)))
             UL=(1-UL)*(10*L-9*LL)+UL*(10*M-9*R);
          else 
             if((Phi_L>=(3./74)) && (Phi_L<(5./6)))
                UL=(1-UL)*((3.*M+6.*L-LL)/8.) + UL*((3.*L+6.*M-R)/8.);
                //UL=(1-UL)*((3./8)*M+(6./8)*L-(1./8)*LL)+UL*((3./8)*L+(6./8)*M-(1./8)*R);
             else 
                if((Phi_L>=(5./6)) && (Phi_L<=1))
                  UL=(1-UL)*M+UL*L;
                else 
                  UL=(1-UL)*L+UL*M;
	
         Der_Conv += (KR*UR-KL*UL)/h;
     }*/
     else
     {
        // Donor_Cell
        if(KR>0)
        {
           if(KL>0) {Der_Conv +=2*(KM*M - KL*L)/h; Sum_coef += 2*(KL - KM);}
           else {Der_Conv += (KR-KL)*M/h; Sum_coef += KL - KR;}
        }
        else 
        {
           if(KL>0) {Der_Conv += (KR*R - KL*L)/h; Sum_coef += KL - KR;}
           else {Der_Conv += 2*(KR*R -KM*M)/h; Sum_coef += 2* (KM - KR);}
        }
        /*// Donor_Cell
        if(KR>0)
        {
           if(KL>0) {Der_Conv +=(KM*M - KLL*L)/h; Sum_coef += KLL - KM;}
           else {Der_Conv += (KR-KL)*M/h; Sum_coef += KL - KR;}
        }
        else 
        {
           if(KL>0) {Der_Conv += (KR*R - KL*L)/h; Sum_coef += KL - KR;}
           else {Der_Conv += (KRR*R -KM*M)/h; Sum_coef += KM - KRR;}
        }*/
	/*// Donor_Cell
        if(KR>0)
        {
           if(KL>0) {Der_Conv +=(KR*M - KL*L)/h; Sum_coef += KL - KR;}
           else {Der_Conv += (KR-KL)*M/h; Sum_coef += KL - KR;}
        }
        else 
        {
           if(KL>0) {Der_Conv += (KR*R - KL*L)/h; Sum_coef += KL - KR;}
           else {Der_Conv += (KR*R -KL*M)/h; Sum_coef += KL - KR;}
        }*/
        /* Substitui texto original apresentado abaixo
        KRABS=fabs(KR);       
        KLABS=fabs(KL);
         FGH[index]=(KR*(R+M) -KL*(L+M))/(hx2); // central differences
        // alpha=0 => central differences,  alpha=1 => full upwind
        FGH[index]=(1-alpha)*du2dx + alpha*( (KR-KRABS)*R + ( (KR+KRABS)-(KL-KLABS) )*M - (KL+KLABS)*L )/(hx2);
        */ 
    } 
}                     

//time handling for auxiliary-velocities

//Compute right hand side for poisson equation
void Navier::CompRHS() 
{ 
   cout<<"CompRHS_begins \n";
   double RHSmax = 0.0;
   int index_max = dj+dk+1;;
   double gridRHSmax[6];
   for(int itfor=0; itfor<List_Fluid_end_zero;itfor++)
   { 
       int index=(*List_Fluid_zero)[itfor];
       //cout<<"index ="<<index<<"\n";
       //cout<<"flag[index] ="<<flag[index]<<"\n";
       //cout<<"List_Fluid_end_zero ="<<List_Fluid_end_zero<<"\n";
       
       if((flag[index] >= FLUID) && (flag[index] < INFLOW) )             
       {
          RHS[index]=  (gridUVW[0][index+1] -gridUVW[0][index] +
	                gridUVW[1][index+dj] -gridUVW[1][index] +
		        gridUVW[2][index+dk] -gridUVW[2][index])/h;
	  //P[index] = 0.0;
	  if(update_pressure) {} //P[index] = 0.0;
	  else
	  {
	     if(fabs(RHS[index]) > RHSmax)
	     { 
	         RHSmax = fabs(RHS[index]);
	         index_max =index;
	         gridRHSmax[0]= gridUVW[0][index];
	         gridRHSmax[1]= gridUVW[1][index];
	         gridRHSmax[2]= gridUVW[2][index];
	         gridRHSmax[3]= gridUVW[0][index+1];
	         gridRHSmax[4]= gridUVW[1][index+dj];
	         gridRHSmax[5]= gridUVW[2][index+dk];
	      }
	      
	  /*if(fabs(RHS[index])> .0001)
	  {
	     int x,y,z;
	     GIJK(index,x,y,z);
	     cout<<x<<","<<y<<","<<z<<"\n";
	     cout<<"RHS ="<<index<<","<<RHS[index]<<","<< gridUVW[0][index+1]<<","<<gridUVW[0][index] <<","<<
	                gridUVW[1][index+dj] <<","<<gridUVW[1][index] <<","<<
		        gridUVW[2][index+dk] <<","<<gridUVW[2][index]<<"\n";
	     int a;
	     cin>>a;
	  }*/		
	  //cout<<" CompRHS -inner if";       
	  }
	}
	else
	{
	 
	 (List_Fluid_end_zero)--;
	 cout<<"List_Fluid_end_zero RHS="<<List_Fluid_end_zero<<","<<itfor<<"\n";
	 cout<<"index ="<<index<<"\n";
         cout<<"flag[index] ="<<flag[index]<<"\n";
	 (*List_Fluid_zero)[itfor]=(*List_Fluid_zero)[List_Fluid_end_zero];
	 //cout<<"List_Fluid_zero[itfor]="<< (*List_Fluid_zero)[itfor]<<"\n";
	 itfor--;
	}  
      
    }
    if(!update_pressure)
    {
       cout<<"RHSmax ="<<RHSmax<<"\n";  
       for(int ifor = 0; ifor < 6; ifor++)  cout<<"gridRHSmax ="<<gridRHSmax[ifor]<<"\n";
       int x,y,z;
       GIJK(index_max,x,y,z);
       cout<<"x,y,z="<<index_max<<","<<x<<","<<y<<","<<z<<"\n";
       cout<<"flags="<<flag[index_max-1]<<","<<flag[index_max+1]<<","<<flag[index_max-dj]<<","<<flag[index_max+dj]<<","<<flag[index_max-dk]<<","<<flag[index_max+dk]<<"\n"; 
       cout<<"pressure="<<P[index_max]<<","<<P[index_max-1]<<","<<P[index_max+1]<<","<<P[index_max-dj]<<","<<P[index_max+dj]<<","<<P[index_max-dk]<<","<<P[index_max+dk]<<"\n";
       double lap_temp = ((P[index_max-1]+P[index_max+1]+P[index_max-dj]+P[index_max+dj]+P[index_max-dk]+P[index_max+dk])- 6* P[index_max])/(h*h);
       cout<<"lap_temp="<<lap_temp<<"\n";
     }
     cout<<"Out of CompRHS \n";
}

//calc new velocities by pressure-gradient correction  
//c�lcula a velocidade final, ap�s o c�lculo da press�o que satisfaz a condi��o de diverg�ncia livre
void Navier::AdapUVW_Time_step() 
{
   
   static double h2 = h * h;
   static double tfconvx2 = 2 * tfconv;
   //static double lim_Viscosity = tfdiff*re*h2/6.0; // Viscosity S->tfdiff = 0.3; S->froude = 1.0;
   double lim_Viscosity =1.0;
   static double constant_part[3] = {4*fabs(g[0])/(froude* h), 4*fabs(g[1])/(froude* h), 4*fabs(g[2])/(froude* h)};
   double deltneu=deltmax ; //S->deltmax � definido no arquivo .nav neste caso = 0.01
   double uvwmax = 0.0;
   int index;
   //cout<<"AdapUVW_List_Fluid_end_zero ="<< List_Fluid_end_zero<<"\n";
   /*for(int itfor=0; itfor<List_Fluid_end_zero;itfor++)
   {
       index=(*List_Fluid_zero)[itfor];
       for(int jt=0;jt<=2;jt++)
       {
           cout<<itfor<<","<<index<<","<<gridUVW[jt][index]<<","<<flag[index]<<"\n";
	   if(abs(gridUVW[jt][index]) > 10.0) exit(1);
       }
   }*/ 
   double deltaPmax = 0.0; 
    for(int itfor=0; itfor<List_Fluid_end_zero;itfor++)
   { 
       index=(*List_Fluid_zero)[itfor];
       double P_cell = P[index];
       for(int jt=0;jt<=2;jt++)
       {
	  int dtemp =d123[jt];
          if(flag[index-dtemp] == FLUID)
             {gridUVW[jt][index]-=/*S->delt* */ (P_cell-P[index-dtemp])/h;  if(fabs((P_cell-P[index-dtemp])/h)>deltaPmax) deltaPmax = fabs((P_cell-P[index-dtemp])/h);}
          else
             if(flag[index-dtemp] == AIR)
                    {gridUVW[jt][index] -=/*S->delt* */ P_cell/h; if(fabs(P_cell/h)>deltaPmax) deltaPmax = fabs(P_cell/h);}
	  double UVW = gridUVW[jt][index];
	  
	  
          deltneu=min(deltneu, tfconvx2/(fabs(UVW)/ h + sqrt( (UVW*UVW)/(h2)+ constant_part[jt]) +1E-15));
	  if(fabs(UVW) > uvwmax)  uvwmax = fabs(UVW);
          if(flag[index+dtemp] == AIR)
	  {
                  gridUVW[jt][index+dtemp] +=/*S->delt* */ P_cell/h;
		  UVW = gridUVW[jt][index+dtemp];
		 /* if(jt == 1 && fabs(UVW)> 0.01)
		     {
		       cout<< "Pression Border - griduvw ="<<index+dtemp<<","<<UVW<<"\n";
		       int a;
		       cin>>a;
		      }
		   if(jt == 0 && fabs(UVW - 0.3)> 0.01)
		     {
		       cout<< "Pression Border - griduvw - 0="<<index+dtemp<<","<<UVW<<"\n";
		       int a;
		       cin>>a;
		      }*/
		  //cout<<itfor<<","<<index<<","<<gridUVW[jt][index+dtemp]<<","<<P_cell<<"\n";
	          //if(abs(gridUVW[jt][index+dtemp]) > 10.0) exit(1);
                  deltneu=min(deltneu, tfconvx2/(fabs(UVW)/ h + sqrt( (UVW*UVW)/(h2)+ constant_part[jt]) +1E-15));
	          if(fabs(UVW) > uvwmax)  uvwmax = fabs(UVW);
	  }
       } 
    }
    uvwmax = fmax(uvwmax_initial,uvwmax);
    double lim_Advection = h/(S->advection_control * uvwmax);
    //cout<<"lim_Advection,h="<<lim_Advection<<","<<uvwmax<<","<<lim_Viscosity<<","<<deltneu<<"\n";
    deltneu=fmin(fmin(deltneu,lim_Advection),lim_Viscosity); 
    if(deltneu<1E-15) 
    {
       printf("Maybe ERROR in parameters: Time-Step = 0.0 ! \n");
       exit(1);    
    }
    S->delt = deltneu;
    *Total_time +=deltneu;
    dt = deltneu;
    cout<<"lim_advection ="<<lim_Advection<<"\n"; 
    cout<<"dt="<< dt <<"\n";
    cout<<"uvwmax="<< uvwmax <<"\n";
    cout<<"deltaPmax="<< deltaPmax <<"\n";
}     
void Navier::Call_Update_Slip_Outflow_Cells() 
{
   
   cout<<"Call_Update_Slip_Outflow_Cells - starting point \n";
   int *neigh_number_pointer, *neigh_pointer;
   int itend =(int)List_Border_Slip_Outflow_Air.size();
   cout<<"itend ="<<itend<<"\n";
   for(int itfor=0; itfor<itend; itfor++) 
   {
            index_cell = List_Border_Slip_Outflow_Air[itfor];
	    
	    if(flag[index_cell]%incrementby2 >= incrementby4)
	    {
	       neighbor_position = flag[index_cell]%incrementby4;
	       index_neighbor = index_cell + cell_to_neighbor_delta[ neighbor_position];
	       if(flag[index_neighbor] >= FLUID)
	       {
		     
		   Update_One_Neighbor_Slip_Outflow_Cells((flag[index_cell] >= OUTFLOW_AIR)  && (flag[index_cell] < AIR_FLUID_LIMIT));
	       }
	       /*if( flag[index_cell] >= AIR_FLUID_LIMIT)
	       {
	           int x,y,z;
	           GIJK(index_cell,x,y,z);
	           cout<<"x,y,z="<<index_cell<<","<<x<<","<<y<<","<<z<<"\n";
	           cout<<"flag[index_cell]-AIR"<<flag[index_cell]<<"\n";  
	        }*/
		//cout<<"Call_Update_Slip_Outflow_Cells - An Air cell processed \n";
	    }
	    else
	    {
	       for(int ifor=0; ifor<6;ifor++) 
	       {
	         grid_temp[ifor]=0.0;
	       }
	       neighbor_code = flag[index_cell]%incrementby4;
	       neigh_number_pointer = &(neighbors_table[neighbor_code][0]);
	       neigh_pointer = neigh_number_pointer + 1;
	       fluid_neighbors = 0;
	       for(int ifor= 0; ifor< *neigh_number_pointer; ifor++)
	       {
	          index_neighbor = index_cell + cell_to_neighbor_delta[*neigh_pointer];
		  if(flag[index_neighbor] >= FLUID)
	          {  
		      Update_Several_Neighbor_Slip_Outflow_Cells((flag[index_cell] >= OUTFLOW_AIR) && (flag[index_cell] < AIR_FLUID_LIMIT));
	          }
		  neigh_pointer++;
		  fluid_neighbors++;
	       }
	       Average_several_fluid_neighbors_contribution();
	    }
	    
   }
  
  // cout<<"Call_Update_Slip_Outflow_Cells - Air cells processed \n";
   itend =(int)List_Slip_Outflow_Fluid.size();
  cout<<"itend ="<<itend<<"\n";
   for(int itfor=0; itfor<itend; itfor++) 
   {
            index_cell = List_Slip_Outflow_Fluid[itfor];
           // if(flag[index_cell] < FLUID)
           // { 
	        if(flag[index_cell]%incrementby2 >= incrementby4)
	        {
		   neighbor_position = flag[index_cell]%incrementby4;
	           index_neighbor = index_cell + cell_to_neighbor_delta[neighbor_position];
		   if(flag[index_neighbor] >= FLUID)
		   { 
		      Update_One_Neighbor_Slip_Outflow_Cells((flag[index_cell] < SLIP_FLUID)  /*&& (flag[index_cell] > AIR_FLUID_LIMIT)*/);
		   }
		   /*if( flag[index_cell] < AIR_FLUID_LIMIT)
	           {
	              int x,y,z;
	              GIJK(index_cell,x,y,z);
	              cout<<"x,y,z="<<index_cell<<","<<x<<","<<y<<","<<z<<"\n";
	              cout<<"flag[index_cell]-FLUID"<<flag[index_cell]<<"\n";
	           }*/
		   //cout<<"Call_Update_Slip_Outflow_Cells - A Fluid cell processed \n";
	        }
	        else
	        {  
		   for(int ifor=0;ifor<6;ifor++) 
	           {
	              grid_temp[ifor]=0.0;
	           }
		   neighbor_code = flag[index_cell]%incrementby4;
	           neigh_number_pointer = &(neighbors_table[neighbor_code][0]);
	           neigh_pointer = neigh_number_pointer + 1;
		   fluid_neighbors = 0;
	           for(int ifor= 0; ifor< * neigh_number_pointer; ifor++)
	           {
	              index_neighbor = index_cell + cell_to_neighbor_delta[*neigh_pointer];
		      if(flag[index_neighbor] >= FLUID)
		      {
		         Update_Several_Neighbor_Slip_Outflow_Cells((flag[index_cell] < SLIP_FLUID) /*&& (flag[index_cell] > AIR_FLUID_LIMIT)*/);
		      }
		      neigh_pointer++;
		      fluid_neighbors++;
	           }
		   Average_several_fluid_neighbors_contribution();
	         }
    }
    //cout<<"Call_Update_Slip_Outflow_Cells - ending point \n";
 }
void Navier:: Update_One_Neighbor_Slip_Outflow_Cells(bool Outflow_or_Slip)
{ 
     
   // cout<<"Update_One_Neighbor_Slip_Outflow_Cells - starting point \n";
     int curr_contact_direction = contact_direction[neighbor_position];
    /* int x,y,z;
     GIJK(index_cell,x,y,z);
     cout<<"x,y,z="<<index_cell<<","<<x<<","<<y<<","<<z<<"\n";
     cout<<"neighbor_position,index_neighbor,curr_contact_direction="<<neighbor_position<<","<<index_neighbor<<","<<curr_contact_direction<<"\n";
     cout<<"index,cell_opposite_face_delta,cell_contact_face_delta,d123="<<index_cell<<","<<cell_opposite_face_delta[neighbor_position]<<","<<cell_contact_face_delta[neighbor_position]
         <<","<<d123[curr_contact_direction]<<"\n";*/
     if(Outflow_or_Slip)
     {
           //cout<<"Label_Outflow ="<<index_cell<<","<<flag[index_cell]<<"\n";int a;cin>>a;
	  double gridUVW_temp = gridUVW[curr_contact_direction][index_cell + cell_contact_face_delta[neighbor_position]];
	  gridUVW[curr_contact_direction][index_cell+ cell_opposite_face_delta[neighbor_position]]= gridUVW_temp;
	  *Total_outflow += (gridUVW_temp*dt);
      }
      else
      {
	   gridUVW[curr_contact_direction][index_cell]= 0.0;
	   gridUVW[curr_contact_direction][index_cell + d123[curr_contact_direction]]=0.0;
      }
      for(int ifor = (curr_contact_direction+1)%3; ifor != curr_contact_direction; ifor =(ifor+1)%3)
      {
           if(flag[index_cell - d123[ifor]] != INFLOW) gridUVW[ifor][index_cell] = gridUVW[ifor][index_neighbor];
           if(flag[index_cell + d123[ifor]] != INFLOW) gridUVW[ifor][index_cell + d123[ifor]] = gridUVW[ifor][index_neighbor + d123[ifor]];
      }
      /*for(int ifor=0;ifor<3;ifor ++)
      {
        cout<<"gridUVW[ifor][index_cell]="<< ifor<<","<<index_cell<<","<<x<<","<<y<<","<<z<<","<<gridUVW[ifor][index_cell]<<","<<gridUVW[ifor][index_cell+d123[ifor]]<<"\n";
      }*/
     // cout<<"Update_One_Neighbor_Slip_Outflow_Cells - ending point \n";
 } 
 void Navier:: Update_Several_Neighbor_Slip_Outflow_Cells(bool Outflow_or_Slip)
{ 
     int curr_contact_direction = contact_direction[neighbor_position];
    // cout<<"Update_Several_Neighbor_Slip_Outflow_Cells - starting point \n";
     if(Outflow_or_Slip)
     {
          cout<<"Label_Outflow ="<<index_cell<<"\n";
          double gridUVW_temp = gridUVW[curr_contact_direction][index_cell + cell_contact_face_delta[neighbor_position]];
	  grid_temp[2*curr_contact_direction + cell_opposite_face_delta_0_1[neighbor_position]] =  gridUVW_temp;
	  *Total_outflow += (gridUVW_temp*dt);
      }
      else
      {
	   gridUVW[curr_contact_direction][index_cell]= 0.0;
	   gridUVW[curr_contact_direction][index_cell + d123[curr_contact_direction]]=0.0;
      }
      int initial_no_contact_face_direction = (curr_contact_direction+1)%3;
      int no_contact_face = 2 * initial_no_contact_face_direction;
      bool *no_contact_pointer = &(not_contact_table[neighbor_code][no_contact_face]);
      for(int ifor = initial_no_contact_face_direction ; ifor != curr_contact_direction; ifor =(ifor+1)%3)
      {
          if(*no_contact_pointer)
	  {
	     grid_temp[no_contact_face] += gridUVW[ifor][index_neighbor];
	  }
	  no_contact_pointer ++;
	  no_contact_face ++;
          if(*no_contact_pointer)
	  {
	     grid_temp[no_contact_face] += gridUVW[ifor][index_neighbor+ d123[ifor]];
	  }
	  if(no_contact_face == 5)
	  {
	    no_contact_pointer -=5;
	    no_contact_face = 0;
	   }
	   else
	   {
	      no_contact_pointer ++;
	      no_contact_face ++;
	   }
      } 
      //cout<<"Update_Several_Neighbor_Slip_Outflow_Cells - ending point \n";
 }
 void Navier::Average_several_fluid_neighbors_contribution()
 { 
       //cout<<"Average_several_fluid_neighbors_contribution-starting point \n";
       int *neigh_number_pointer;
       int *neigh_pointer ;
       neigh_number_pointer = &(neighbors_table[63 - neighbor_code][0]);
       neigh_pointer = neigh_number_pointer + 1;
       //cout<<"neighbors_table[63 - neighbor_code][0] ="<<neighbors_table[63 - neighbor_code][0]<<"\n";
       //cout<<"fluid_neighbors ="<<fluid_neighbors<<"\n";
       //cout<<"index_cell="<<index_cell<<"\n";
       //int x,y,z;
       //GIJK(index_cell,x,y,z);
       //cout<<"x,y,z="<<x<<","<<y<<","<<z<<"\n";
       //cout<<"flag[index_cell]="<<flag[index_cell]<<"\n";
       if(fluid_neighbors == 0) return;
       for(int ifor= 0; ifor< * neigh_number_pointer; ifor++)
       {
	   int face_index = *neigh_pointer;
	   int index_no_contact_face = index_cell + cell_contact_face_delta[face_index];
	   gridUVW[face_index/2][index_no_contact_face]= grid_temp[face_index]/ fluid_neighbors; 
	   neigh_pointer++; 
       }
       //cout<<"Average_several_fluid_neighbors_contribution-ending point \n";
 }     	      
Navier::~Navier() { }
